
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! +                                                           +
! +  glimmer_utils.f90 - part of the Glimmer-CISM ice model   + 
! +                                                           +
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 
! Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010
! Glimmer-CISM contributors - see AUTHORS file for list of contributors
!
! This file is part of Glimmer-CISM.
!
! Glimmer-CISM is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 2 of the License, or (at
! your option) any later version.
!
! Glimmer-CISM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Glimmer-CISM.  If not, see <http://www.gnu.org/licenses/>.
!
! Glimmer-CISM is hosted on BerliOS.de:
! https://developer.berlios.de/projects/glimmer-cism/
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

!> Module containing utility code for GLIMMER.
module glimmer_utils


  use glimmer_global

  implicit none

  interface array_bcs
    module procedure array_bcs1d,array_bcs2d
  end interface

  interface check_conformal
    module procedure check_conformal_2d_real
  end interface

contains

  !> Returns the value of a 1D array location,checking first for the boundaries.
  !!
  !! the location is wrapped around the array boundaries until it falls within the array
  !! \author The value of the location in question.
  real(rk) function array_bcs1d(array,i)

    ! Arguments

    real(rk),dimension(:),intent(in) :: array !< The array to be indexed.
    integer,intent(in)               :: i     !< The location to be extracted.

    ! Internal variables

    integer :: n,ii

    n=size(array)
    ii=i

    if ((i<=n).and.(i>=1)) then
      array_bcs1d=array(i)
    endif

    do while (ii>n)
      ii=ii-n
    enddo

    do while (ii<1)
      ii=ii+n
    enddo

    array_bcs1d=array(ii)

  end function array_bcs1d

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Returns the value of a 1D array location,checking first for the boundaries.
  !!
  !! the location is wrapped around the array boundaries until it falls within the array
  !! as array_bcs1d but for polar boundary conditions
  !! \author The value of the location in question.
  real(rk) function array_bcs_lats(array,i)


    ! Arguments

    real(rk),dimension(:),intent(in) :: array !< The array to be indexed.
    integer,intent(in) :: i !< The location to be extracted.

    ! Internal variables

    integer :: n,ii

    n=size(array)
    ii=i

    if ((i<=n).and.(i>=1)) then
      array_bcs_lats=array(i)
      return
    endif

    if (ii>n) then
      ii=2*n-ii
      array_bcs_lats=-180.0+array(ii)
    endif

    if (ii<1) then
      ii=1-ii
      array_bcs_lats=180.0-array(ii)
    endif

  end function array_bcs_lats

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Returns the value of an array 
  !! location, checking first for the boundaries. 
  !! Over-the-pole boundary conditions are implemented here.
  !! \return The value of the location specified.
  real(rk) function array_bcs2d(array,i,j)

    ! Arguments

    real(rk),dimension(:,:),intent(in) :: array !< Array to be indexed
    integer,intent(in) :: i !< The location to be extracted    
    integer,intent(in) :: j !< The location to be extracted

    ! Internal variables

    integer :: nx,ny,ii,jj

    nx=size(array,1) ; ny=size(array,2)

    if ((i>=1).and.(i<=nx).and.(j>=1).and.(j<=ny)) then
      array_bcs2d=array(i,j)
      return
    endif

    ii=i ; jj=j

    if (jj>ny) then
      jj=2*ny-jj
      ii=ii+nx/2
    endif

    if (jj<1) then
      jj=1-jj
      ii=ii+nx/2
    endif

    do while (ii>nx) 
      ii=ii-nx
    enddo

    do while (ii<1)
      ii=ii+nx
    enddo  

    array_bcs2d=array(ii,jj)

  end function array_bcs2d

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Adjusts array location indicies
  !! so that they fall within the domain.
  subroutine fix_bcs2d(i,j,nx,ny)


    integer,intent(inout) :: i !< The location of interest
    integer,intent(inout) :: j !< The location of interest
    integer,intent(in) :: nx  !< The size of the domain (number of points in each direction)
    integer,intent(in) :: ny  !< The size of the domain (number of points in each direction)

    if ((i>=1).and.(i<=nx).and.(j>=1).and.(j<=ny)) return

    if (j>ny) then
      j=2*ny-j
      i=i+nx/2
    endif

    if (j<1) then
      j=1-j
      i=i+nx/2
    endif

    do while (i>nx) 
      i=i-nx
    enddo

    do while (i<1)
      i=i+nx
    enddo  

  end subroutine fix_bcs2d

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Checks that two arrays are of the same size.
  subroutine check_conformal_2d_real(array1,array2,label)

    use glimmer_log


    real(rk),dimension(:,:),intent(in) :: array1 !< The array 1 to be checked
    real(rk),dimension(:,:),intent(in) :: array2 !< The array 2 to be checked
    character(*),intent(in),optional :: label    !< Optional label, to facilitate bug tracking if the check fails.

    if ((size(array1,1)/=size(array2,1)).or.(size(array1,2)/=size(array2,2))) then
      if (present(label)) then
        call write_log('Non-conformal arrays. Label: '//label,GM_FATAL,__FILE__,__LINE__)
      else
        call write_log('ERROR: Non-conformal arrays. No label',GM_FATAL,__FILE__,__LINE__)
      endif
    endif

  end subroutine check_conformal_2d_real

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> compute horizontal sum for each vertical level
  !!
  !! Calculates the sum of a given three-dimensional field at each
  !! level. The vertical coordinate of the input is the first index of
  !! the array.
  !! \return
  !! A one-dimensional array of the same size as the first dimension of
  !! inp is returned, containing the sum of inp for 
  !! each level.
  function hsum(inp)


    implicit none

    real(dp),dimension(:,:,:),intent(in) :: inp !< The input array. The first index is the vertical, the othe two horizontal.
    real(dp),dimension(size(inp,dim=1))  :: hsum
  
    integer up

    do up=1,size(inp,dim=1)
       hsum(up) = sum(inp(up,:,:))
    end do

  end function hsum

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> compute horizontal sum of a 2x2 horizontal mesh for each vertical level
  !!
  !! Calculates the sum of a given three-dimensional field at each
  !! level. The vertical coordinate of the input is the first index of
  !! the array.
  !! \return
  !! A one-dimensional array of the same size as the first dimension of
  !! inp is returned, containing the sum of inp for 
  !! each level.  
  function hsum4(inp)

    implicit none

    real(dp),dimension(:,:,:),intent(in) :: inp !< The input array. The first index is the vertical, the other two horizontal.
    real(dp),dimension(size(inp,dim=1))  :: hsum4
  
    hsum4(:) = inp(:,1,1) + inp(:,2,1) + inp(:,1,2) + inp(:,2,2)

  end function hsum4

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Calculates the sum of a given two-dimensional field along one axis.
  !! Within GLIMMER, this function calculates the mean vertical profile
  !! in a 2D vertical slice. 
  !! \return
  !! A one-dimensional array of the same size as the first dimension of
  !! inp is returned, containing the sum of inp for 
  !! each row.
  function lsum(inp)


    implicit none

    real(dp),dimension(:,:), intent(in) :: inp !< Input array
    real(dp),dimension(size(inp,dim=1)) :: lsum
    
    lsum = sum(inp(:,:),dim=2)

  end function lsum

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Tridiagonal solver. All input/output arrays should have the 
  !! same number of elements.
  subroutine tridiag(a,b,c,x,y)


    real(dp),dimension(:) :: a !< Lower diagonal; a(1) is ignored.
    real(dp),dimension(:) :: b !< Centre diagonal
    real(dp),dimension(:) :: c !< Upper diagonal; c(n) is ignored.
    real(dp),dimension(:) :: x !< Unknown vector
    real(dp),dimension(:) :: y !< Right-hand side

    real(dp),dimension(size(a)) :: aa
    real(dp),dimension(size(a)) :: bb

    integer :: n,i

    n=size(a)

    aa(1) = c(1)/b(1)
    bb(1) = y(1)/b(1)

    do i=2,n
       aa(i) = c(i)/(b(i)-a(i)*aa(i-1))
       bb(i) = (y(i)-a(i)*bb(i-1))/(b(i)-a(i)*aa(i-1))
    end do
    
    x(n) = bb(n)

    do i=n-1,1,-1
       x(i) = bb(i)-aa(i)*x(i+1)
    end do

  end subroutine tridiag

end module glimmer_utils
